Br = 1.465;
nu0 = 4*pi*10^(-7);
effector.radius = 15e-3;
effector.diameter = 2*effector.radius;
effector.height = 30e-3;
effector.volume = pi*(effector.radius^2)*effector.height;
effector.magDirection = [0;0;1];
effector.postion.Cur = [0;0;0];
effector.attitude0 = [pi/2;0;0];

% driver.radius = 15e-3;
% driver.diameter = 2*driver.radius;
% driver.height = 30e-3;
% driver.volume = pi*(driver.radius^2)*driver.height;
% driver.magDirection = [0;0;1];
% driver.position.Cur = [0;100e-3;0];
% driver.attitude0 = [pi/2;0;0];
diffBmap = [];
for z = linspace(40e-3,60e-3,20)
    for y = linspace(-30e-3,30e-3,61)
        position = [0;y;z];
        diffBmap = [diffBmap [position;diffMapFun(effector,position,'dip')*[0;0;1]]];
    end
end
z = linspace(25e-3,50e-3,20);
y = linspace(-30e-3,30e-3,61);
[xx,yy]=meshgrid(z,y);
diff_x = diffBmap(4,:);
diff_y = diffBmap(5,:);
diff_z = diffBmap(6,:);
diff_y = reshape(diff_y,61,20);
%% 
figure(1)
plot3(xx,yy,diff_y)
figure(2)
set(gca,'FontSize',60,'Fontname', 'Times New Roman');
surfc(1000*xx,1000*yy,diff_y)
xlabel('z(mm)');
ylabel('y(mm)')
zlabel('\partial B_z/\partial y(T/m)');
set(gca,'FontSize',40,'Fontname', 'Times New Roman');
hold on
contour3(1000*xx,1000*yy,diff_y,61) 
ContourZLevel = -3;
%%
function diff = diffMapFun(effector,position,method)
% 计算梯度矩阵，其中efferctor是磁源头，method是选用的方法
Br = 1.465;
dx = 0.01e-3;dy =0.01e-3;dz=0.01e-3;
driver.position.Cur = position;
driver.position.xPre  =  driver.position.Cur - [dx/2;0;0];
driver.position.xNext =  driver.position.Cur + [dx/2;0;0];
driver.position.yPre  =  driver.position.Cur - [0;dy/2;0];
driver.position.yNext =  driver.position.Cur + [0;dy/2;0];
driver.position.zPre  =  driver.position.Cur - [0;0;dz/2];
driver.position.zNext =  driver.position.Cur + [0;0;dz/2];
switch method
    case 'ana'
        effector.magInduction.Ana.xPre  = PMCir3D(Br,effector.diameter,effector.height,driver.position.xPre);
        effector.magInduction.Ana.xNext = PMCir3D(Br,effector.diameter,effector.height,driver.position.xNext);
        effector.magInduction.Ana.yPre  = PMCir3D(Br,effector.diameter,effector.height,driver.position.yPre);
        effector.magInduction.Ana.yNext = PMCir3D(Br,effector.diameter,effector.height,driver.position.yNext);
        effector.magInduction.Ana.zPre  = PMCir3D(Br,effector.diameter,effector.height,driver.position.zPre);
        effector.magInduction.Ana.zNext = PMCir3D(Br,effector.diameter,effector.height,driver.position.zNext);
        effector.magInduction.Ana.diffMatrix = [(effector.magInduction.Ana.xNext-effector.magInduction.Ana.xPre)';...
            (effector.magInduction.Ana.yNext-effector.magInduction.Ana.yPre)';...
            (effector.magInduction.Ana.zNext-effector.magInduction.Ana.zPre)'  ];
        diff = effector.magInduction.Ana.diffMatrix/dx;
    case 'dip'
        effector.magInduction.Dip.xPre  = PMDipCir(effector.radius,effector.height,driver.position.xPre ,'z');
        effector.magInduction.Dip.xNext = PMDipCir(effector.radius,effector.height,driver.position.xNext,'z');
        effector.magInduction.Dip.yPre  = PMDipCir(effector.radius,effector.height,driver.position.yPre ,'z');
        effector.magInduction.Dip.yNext = PMDipCir(effector.radius,effector.height,driver.position.yNext,'z');
        effector.magInduction.Dip.zPre  = PMDipCir(effector.radius,effector.height,driver.position.zPre ,'z');
        effector.magInduction.Dip.zNext = PMDipCir(effector.radius,effector.height,driver.position.zNext,'z');
        effector.magInduction.Dip.diffMatrix = [(effector.magInduction.Dip.xNext-effector.magInduction.Dip.xPre)';...
            (effector.magInduction.Dip.yNext-effector.magInduction.Dip.yPre)';...
            (effector.magInduction.Dip.zNext-effector.magInduction.Dip.zPre)'   ];
        diff = effector.magInduction.Dip.diffMatrix/dx;
end
end
